Modelagem geoestatística da profundidade do solo em área de produção silvicultural
Horst, Taciara Zborowski
1 de novembro de 2017
1 Objetivo
Realizar a modelagem geoestatística da profundidade do solum área de produção silvicultural
2 Caracterização da área de estudo
Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de Pinus taeda L. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.
A área possui Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).
Figure 2.1: Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).
3 Amostragem
A coleta de dados foi realizada em 94 pontos (Figura 1 b) alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função clhs do pacote clhs para R. Foram consideradas como variáveis condicionantes da amostragem a elevação, profundidade do vale, índice de umidade topográfico, nivel base da rede de drenagem e declividade as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).
Em cada ponto amostral foi medida a profundidade do solum (PF). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras.
Cada ponto de amostragem representa uma área de aproximadamente 100 m2 no campo. Apesar do suporte amostral ser areal, considerei-os como sendo em ponto.
Os dados foram armazenados no objeto pontos definido como SpatialPolygonsDataFrame com as coordenadas projetadas em WGS84.
pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)
pontos$PFd <- pontos$PF/10
sp::coordinates(pontos) <- c('X' , 'Y')
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sp::proj4string(pontos) <- wgs84utm22sA localização dos pontos no espaço pode ser observada a seguir:
mapview(pontos, zcol = "PFd")4 Modelagem geoestatística da profundidade do Solum
4.1 Dados de profundidade do solum
Como uma análise preliminar avaliei a distribuição de frequências dos dados de profundidade do solum (Figura 2). Os dados apresentam uma distribuição bimodal, com valores entre 1 e 10 e média de 5,86 dm.
Figure 4.1: Histograma de frequências da profundidade do solum
summary(pontos$PFd)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 2.000 6.000 5.858 10.000 10.0004.2 Modelo linear misto de variação espacial
Para modelagem geoestatística dos dados foi utilizado o modelo linear misto de variação espacial denotado por
\[Y(\boldsymbol{s}_i) = Z(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i) = \boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta} + B(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i)\]
Para utilizar esse modelo foi necessário supor que os dados são uma realização de um campo aleatório \(Y(\boldsymbol{s}_i)\) com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos estocásticos e erro aleatório independente.
\(Z(\boldsymbol{s}_i)\) ou sinal possui dois componentes. O primeiro (efeito fixo) \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) representa os efeitos de origem desterminística, que relaciona a dependência entre a variável e as covariáveis.
O segundo componente do sinal (efeito aleatório), \(B(\boldsymbol{s}_i)\), um campo aleatório Gaussiano estacionário não-observável, descrito por sua função de média e função de covariância.
\(\varepsilon(\boldsymbol{s}_i)\) é o erro (ou ruído), descrito por uma distribuição Gaussina de probabilidade, cujo parâmetro desconhecido de escala é \(\tau\).
4.2.1 Efeito fixo
Para selecionar os efeitos fixos do modelo, verifiquei a relação linear entre as variáveis topográficas e de produtividade e a profundidade do solum.
4.2.1.1 Covariáveis topográficas
As covariáveis topográfricas foram derivadas do modelo digital de elevação disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta reamostragem no software SAGA GIS.
Os planos de informação utilizados foram importados pela ferramenta raster::rastere armazenados em .tif nos objetos RasterLayer (Figura 3).
DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
VD <- raster::raster("../data/Covars/VD.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
Figure 4.2: Covariáveis topográficas utilizadas como preditoras no modelo floresta aleatória
A partir da função raster::extract extraí os valores de cada objeto na localização de cada observação contida no objeto espacial pontos.
pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$VD <- raster::extract(VD, pontos)
pontos$TWI <- raster::extract(TWI, pontos)4.2.1.2 Covariável de produção
Além dos atributos de terreno, foi utilizado um índice de produtividade disponibilizado pela empresa.
Esses índices são utilizados para ordenamento da produção e provém da amostragem de parcelas fixas. A área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local. Cada PIC é classificada em função da média de altura das 100 árvores de maior perímetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrementos anuais é estabelecido o índices de sítio (Sitio). Esse índice então é então atribuido à todo o talhão.
As parcelas da área são classificadas em 4 níveis, em que 1 corresponde a melhor e 4 a pior produtividade (Figura 4).
Utilizei a função raster::shapefile para carregar o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM) e a função raster::extract para extraír os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial pontos.
Sitio <- raster::raster("../data/Covars/Sitio.tif")
Sitio <- as.factor(Sitio)
sp::proj4string(Sitio) <- wgs84utm22s
pontos$Sitio <- raster::extract(Sitio, pontos) %>% as.factor()
pontos <- na.omit(pontos)ProdutividadePistola <-
raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
sp::spTransform(wgs84utm22s)
ProdutividadePistola$Sitio <- as.factor(ProdutividadePistola$Sitio)
sp::spplot(
ProdutividadePistola, scales = list(draw = T),
main = 'Classificação dos sítios e localização dos pontos') +
lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()
Figure 4.3: Classificação dos sítios e localização dos pontos na área de estudo
4.2.1.3 Avaliação da relação linear entre as variáveis
Considerando que apenas a covariável DECLI apresenta uma correlação significativa com a profundidade do solum, considerei somente a declividade como variável de efeito fixo, não-constantes, na área.
lm(PFd ~ DECLI + ELEV + VD + Sitio, data = pontos@data) %>% summary()##
## Call:
## lm(formula = PFd ~ DECLI + ELEV + VD + Sitio, data = pontos@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.235 -2.569 -0.349 2.922 5.235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.249e+01 4.116e+01 0.303 0.76223
## DECLI -1.517e+01 5.341e+00 -2.840 0.00557 **
## ELEV -4.219e-03 4.458e-02 -0.095 0.92481
## VD -8.306e-04 5.096e-02 -0.016 0.98703
## Sitio2 1.059e-01 1.089e+00 0.097 0.92270
## Sitio3 -2.024e+00 1.105e+00 -1.832 0.07017 .
## Sitio4 -1.257e+00 9.954e-01 -1.263 0.20998
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.381 on 91 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.1666, Adjusted R-squared: 0.1116
## F-statistic: 3.031 on 6 and 91 DF, p-value: 0.0095254.3 Variograma amostral
O variograma amostral (Figura 7) foi computado através da função georob::sample.variogram. O estimador para semivariãncia foi Metheron (método dos momentos). Para a obtenção dos parâmetros utilizei um corte de 66% da distância máxima entre os pontos, armazenada no objeto distmax, excluindo os pares de longo alcance.
distmax <-dist(pontos@coords) %>% max() /3
limites <- seq(0, distmax, length.out = 20)
vario<- georob::sample.variogram(PFd ~ DECLI,
data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = "matheron") %>%
plot(ylab = 'Seminvariância', xlab = 'Distância de separação (m)', annotate.npairs = TRUE, main = "Semivariograma")
Figure 4.4: Variograma amostral
4.4 Ajuste da função que descreve o varigrama amostral
O método de ajuste empregado no variograma amostral foi o de quadrados mínimos não-lineares ponderados, com ajuste de um modelo exponencial, com ponderação definida conforme o método de “Cressie”. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = "BFGS"). A função resultante desse processo foi armazenada no objeto vario_fit.
vario_fit <-
georob::fit.variogram.model(
vario, variogram.model = 'RMexp', param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", method = "BFGS")##
##
## variance snugget nugget scale
## Variogram parameters 7.666473e-122 0.000000e+00 1.458980e-65 9.136816e+14summary(vario_fit)##
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp",
## param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie",
## method = "BFGS")
##
## Convergence in 97 function and 34 Jacobian/gradient evaluations
##
## Residual Sum of Squares: 58.86111
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -4.0586 -1.8355 -0.6017 1.5709 2.8500
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 12.43754 9.10081 17.00
## snugget(fixed) 0.00000 NA NA
## nugget 1.28807 0.06083 27.28
## scale 72.24013 52.67343 99.08O ajuste do modelo exponencial ao variograma amostral é mostrado na Figura 8. A curva ajustada passa próximo ao centro de massa dos vinte pontos do variograma amostral.
plot(vario, xlab = 'Distância de separação (m)', ylab = 'Semivariância', annotate.npairs = TRUE)
lines(vario_fit, col = "red", lty = 'dashed')
Figure 4.5: Variograma amostral (em preto) e função exponencial a ele ajustado (vermelho)
Considerando que durante a obteção dos dados de campo as informações de profundidade do solum os valores foram obtidos em decímetros, considerei que a variância do erro de medida é igual a 0,5 dm devido ao arredondamento dos valores.
Assim, foi possível discretizar a variância não explicada em erro de medida corresponde 0,25 do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget. A função resultante desse processo foi armazenado no objeto vario_fit_error.
nugget <- 0.25
vario_fit_error <- georob::georob(
PFd ~ DECLI, pontos, locations = ~ X + Y, variogram.model = 'RMexp',
param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']],
nugget = nugget,
snugget = vario_fit$variogram.object[[1]]$param[['nugget']] - nugget,
scale = vario_fit$variogram.object[[1]]$param[['scale']]),
fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE),
tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)##
## Call:georob::georob(formula = PFd ~ DECLI, data = pontos, locations = ~X +
## Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]],
## nugget = nugget, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] -
## nugget, scale = vario_fit$variogram.object[[1]]$param[["scale"]]),
## fit.param = georob::default.fit.param(nugget = FALSE, snugget = TRUE),
## tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
##
## Tuning constant: 1000
##
## Convergence in 8 function and 6 Jacobian/gradient evaluations
##
## Estimating equations (gradient)
##
## xi scale
## Gradient : 1.069608e-03 -8.586508e-04
##
## Maximized restricted log-likelihood: -263.634
##
## Predicted latent variable (B):
## Min 1Q Median 3Q Max
## -5.8300 -2.9564 -0.1296 3.6171 5.4093
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -0.1847714 -0.0623897 0.0009239 0.0657684 0.1921513
##
## Standardized residuals:
## Min 1Q Median 3Q Max
## -2.29890 -0.79556 0.01281 0.88223 1.98182
##
##
## Gaussian REML estimates
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 10.696 7.607 15.04
## snugget(fixed) 1.038 NA NA
## nugget(fixed) 0.250 NA NA
## scale 48.614 24.323 97.16
##
##
## Fixed effects coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2931 0.7609 9.585 8.05e-16 ***
## DECLI -12.1389 5.0173 -2.419 0.0174 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sqrt(nugget)): 0.5
##
## Robustness weights:
## All 102 weights are ~= 1.A comparação entre as funções vario_fit_error e vario_fit está na figura 9.
plot(vario)
lines(vario_fit, col = "blue")
lines(vario_fit_error, col = "red")
Figure 4.6: Variograma amostral (em preto) do modelo linear da profundidade do solum e a função exponencial a ele ajustada (azul) e a função ajustada com erro de medida fixo
A função vario_fit_error foi utilizada para a predição espacial. Para isso foi criado um grid de predição abrangendo a área de estudo. A partir da função raster::extract foram extraídos os valores do objeto DECLI e armazenados no novo conjunto de dados grid.
grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
grid$DECLI <- raster::extract(DECLI, grid)
grid<-
sp::SpatialPointsDataFrame(
coords = grid@coords,
data = data.frame(grid),
proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)4.5 Predição espacial da resposta do modelo linear misto
Utilizei a função predict para aplicar a função vario_fit_erroraos pontos de toda área armazenados no objeto grid
pred_ponto <- predict(
vario_fit_error, newdata = grid, type= "response", signif = 0.95,
control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE
spplot(pred_ponto)
Figure 4.7: Mapas de predição - saída extendida do georob
trend é a média espacial dos efeitos deterministicos (declividade), suavisados.
A krigagem tenta minimizar o valor esperado do erro entre o valor predito e o medido. Ainda assim o intervalo de predição ainda é grande
A raiz do SE é erro padrão da krigagem. O quadrado do SE é a variância dos erros de predição.
at <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range(na.rm = TRUE)
at <- seq(at[1], at[2], length.out = 20)
spplot(pred_ponto, zcol="lower", at=at, main="Limite inferior")
spplot(pred_ponto, zcol="pred", at=at, main="Realização mais provável")
spplot(pred_ponto, zcol="upper", at=at, main="Limite superior")
4.5.1 Avaliação da qualidade da predição
Para a validação dos resultados utilizei o método de validação cruzada. Este método é realizado em etapas de validação com a partição do conjunto de dados aleatóriamente em subconjuntos (\(k\)). Em cada etapa de validação um dos conjuntos é utilizado para validar o modelo calibrado com \(k-1\) subconjuntos. Esse procedimento é repetido até que todos subconjuntos \(k\) tenham sido utilizados como conjunto de validação do modelo. A partir das predições realizadas com cada subconjunto \(k\) são calculados os erros para avaliar a qualidade das predições.
Para isso utilizei a função cv::georob. Defini o número o número \(k\) de subconjuntos em que o conjunto de dados pontos foi particionado definindo o parâmetro nset = 101. Assim, em cada etapa da validação 101 pontos foram utilizados para a validação e 1 para a calibração do modelo linear misto.
validacao <- georob::cv(vario_fit_error, nset = 101)
summary(validacao)##
## Statistics of cross-validation prediction errors
## me mede rmse made qne msse medsse
## -0.04291 0.03903 3.29989 4.27377 3.35070 1.08072 0.76194
## crps
## 1.924581 - sum((validacao$pred$data - validacao$pred$pred)^2) / sum((validacao$pred$data - mean(validacao$pred$data))^2)## [1] 0.1323907plot(validacao)
Como os dados utilizados para a validação foram os mesmo utilizados para validar o modelo, e, que as amostras foram obtidas intencionalmente, é provavel que esta avaliação seja otimista.